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SUMMARY 

During  the  period  of  performance  a  program  of  theoretical  studies  was  carried  out  with  the 
objective  of  quantifying  the  requirements  for  controlled  precipitation  of  radiation  belt  particles 
in  order  to  mitigate  enhanced  space  particle  effects,  to  address  critical  issues  of  antenna-plasma 
coupling,  and  to  develop  a  model  of  the  effective  radiated  power  of  both  dipole  and  loop  an¬ 
tennas  throughout  the  magnetosphere.  In  carrying  out  this  effort,  Stanford  University  worked 
in  conjunction  with  AFRL/VSB  scientists.  As  a  result  of  these  studies,  Stanford  University  has 
developed  a  model  which  for  the  first  time  quantifies  the  radiated  electromagnetic  power  in  the 
1-10  kHz  frequency  range  which  is  necessary  to  reduce  the  lifetimes  of  energetic  electrons  in 
the  radiation  belts  by  an  order  of  magnitude.  In  addition,  a  method  was  developed  to  control 
the  ion  and  electron  sheaths  which  surround  dipole  antennas  in  a  plasma.  Without  such  control, 
the  sheath  effects  strongly  diminish  the  radiation  efficiency  of  dipole  antennas,  and  much  larger 
transmitter  power  is  required  to  produce  the  electromagnetic  wave  fields  necessary  to  control 
the  radiaton  belts. 
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FINAL  TECHNICAL  REPORT 


1.  Contract  Purpose 

The  overall  purpose  of  this  work  was  to  quantify  the  requirements  for  controlled  precipitation 
of  radiation  belt  particles  in  order  to  mitigate  enhanced  space  particle  effects,  to  address  critical 
issues  of  antenna-plasma  coupling,  and  to  develop  a  model  of  the  effective  radiated  power  from 
dipole  antennas  throughout  the  magnetosphere.  In  carrying  out  this  effort,  Stanford  University 
worked  closely  with  AFRL/VSB  scientists.  Stanford  University  concentrated  on:  1)  developing 
a  model  to  quantify  the  radiated  electromagnetic  power  in  the  1-10  kHz  frequency  range  which 
is  necessary  to  reduce  the  lifetimes  of  energetic  electrons  in  the  radiation  belts  by  an  order  of 
magnitude,  and  2)  developing  a  method  to  control  the  ion  and  electron  sheaths  which  surround 
dipole  antennas  in  a  plasma, 

2.  Period  of  Performance 

The  period  of  performance  under  this  contract  extended  from  September  30,  1999,  through 
December  31,  2002. 

3.  Work  Provided 

During  the  period  of  performance  of  this  grant,  Stanford  University  carried  out  the  following 
tasks: 

Task  1:  Stanford  developing  a  model  to  quantify  the  radiated  electromagnetic  power  in  the 
1-10  kHz  frequency  range  which  is  necessary  to  reduce  the  lifetimes  of  energetic  electrons  in 
the  radiation  belts  by  an  order  of  magnitude. 

Task  2:  Stanford  developed  a  method  to  control  the  ion  and  electron  sheaths  which  surround 
dipole  antennas  in  a  plasma. 


2 


Task  3:  Stanford  developed  a  method  to  calculate  the  near  field  and  the  radiated  electromag¬ 
netic  fields  from  dipole  antennas  in  the  radiation  belts. 

Task  3:  Stanford  performed  grant  management,  reporting,  and  technical  overview  activities 
of  work. 

6,  Results  of  Research  Effort 

Detailed  results  of  the  Stanford  research  effort  are  reported  in  Appendix  A. 

7,  List  of  personnel  contributing  to  report 

The  list  of  Stanford  University  scientists  and  engineers  who  contributed  to  the  work  reported  in 
this  document  is  as  follows: 

Tim  Bell 
Jacob  Bortnik 
Tim  Chevalier 


Umran  Inan 


APPENDIX  A. 

Progress  Report  for  AFRL  Grant  F49620-99-1-0339  Entitled: 
Controlled  Precipitation  of  Radiation  Belt  Particles 
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AFRL  Grant  F49620-99-1-0339  has  been  carried  out  in  cooperation  with  AFRL/VSB. 

Below  in  Sections  A,  B,  and  C  we  list  the  significant  accomplishments  resulting  from 
work  carried  out  predominantly  at  Stanford  University  STAR  Laboratory.  Joint  work 
with  AFRL/VSB  personnel  is  discussed  in  Section  D. 

A.  Simulation  Models. 

The  first  Stanford  goal  has  been  the  determination  of  the  radiation  efficiency  and  input 
impedance  of  dipole  antennas  in  the  magnetosphere,  including  the  effects  of  the  plasma  sheath 
surrounding  each  antenna  element.  An  antenna  operating  at  ELF/VLF  frequencies  in  the 
magnetosphere  will  radiate  waves  with  wavelengths  of  many  kilometers.  At  the  same  time  the 
sheath  surrounding  this  antenna  will  have  a  scale  which  ranges  from  centimeters  to  meters  as  the 
antenna  voltage  is  varied.  To  completely  resolve  this  space  from  antenna  to  far  field  using  a 
fixed  cell  size  of  1  cm  in  a  FDTD  simulation  would  require  1000003  cells.  This  requirement  is 
far  greater  than  that  which  even  the  most  powerful  of  contemporary  supercomputers  can  satisfy. 
Thus  it  is  necessary  to  partition  the  simulation  space  into  regions  in  which  different  scale  sizes 
prevail.  As  illustrated  in  Figure  1,  we  have  chosen  to  partition  the  simulation  space  into  three 
separate  regions:  1)  the  plasma  sheath  region,  extending  out  to  a  few  meters  from  the  antenna,  in 
which  either  electrons  or  ions  predominate  according  to  the  sign  of  the  antenna  voltage,  and  in 
which  these  particles  are  accelerated  to  high  velocities  by  the  quasi-electrostatic  fields;  2)  the 
nonlinear  warm  plasma  region,  extending  from  a  few  meters  out  to  100  m  from  the  antenna,  in 
which  the  plasma  is  essentially  neutral  but  the  acceleration  due  to  the  electric  and  magnetic  fields 
produced  by  the  antenna  currents  is  large  enough  to  preclude  the  use  of  a  linear  cold  plasma 
formulation;  and  3)  the  cold  plasma  region  in  which  the  plasma  is  essentially  neutral  and  particle 
energization  is  small. 

In  order  to  make  progress  on  this  topic  we  have  developed  a  series  of  simulation  models, 
each  of  which  applies  to  one  of  the  regions  shown  in  Figure  1 .  To  date  we  have  tested 
and  verified  each  model.  Our  ultimate  goal  is  to  combine  the  simulation  models  to 
produce  a  complete  description  of  the  particles  and  fields  within  the  entire  space.  The 
most  complicated  simulation  model  is  that  which  applies  to  region  2,  the  hot  plasma 
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region.  The  simulation  here  involves  a  3-D  fluid  model  of  the  dynamics  of  electrons  and  ions 
under  the  influence  of  the  electromagnetic  and  quasi-electrostatic  fields  produced  by  the  currents 
and  charges  on  the  dipole  antennas.  We  discuss  this  model  first  below.  The  other  models  are 
described  thereafter. 


Figure  1,  Regions  of  Simulation  Space 


A.l  Full  EM  Fluid  Code 

The  Full  Electromagnetic  Fluid  Code  is  applied  to  region  2  of  Figure  1  and  includes  Maxwell’s 
equations  for  modeling  the  electromagnetic  and  electrostatic  portion  of  the  code,  and  a  fluid 
model  for  determining  the  first  three  moments  of  the  Boltzmann  equation  including  conservation 
of  number  density,  momentum,  and  energy  for  one  species  of  ion  and  electrons.  The  equations 
are  listed  as  follows: 


3 


— + V-(pu)=  Sp 
dt 

■^p-+V-(puu)  =  — V*n  +  gr|(E+ux  2?)+uSp  +Sn 


Q_ 

dt 


1  2  1 
-pu  + - —p 

2  y  -1 


Vx  E  =  - 


m 

dt 


dt 


+  ^-(^‘Pu2u+~Yi?u+ii-n  +  Z)=  <grqu-E+^-w2Sp  +u-Sm  +  Sr 


3=T)eqeUe  +T^U, 


Sa  =  sources 

p  =  mass  density 

p  =  pressure 

7]  =  number  Density 

u  =  Bulk  Velocity  of  Particles 

E  =  electric  Field 

B  =  Bq  +  B  =  Static  Magnetic  Field  +  Wave  Magnetic  Field 
m  =  mass 

n  =  pressure  tensor 
L  =  heat  flux 
J  =  current  density 

A  number  of  simplifications  can  be  made  depending  upon  where  in  the  magnetosphere  the 
antenna  will  be  operating.  Since  it  is  envisioned  that  the  antenna  will  be  operating  on  a 
spacecraft  near  the  magnetic  equatorial  plane  at  L=2  or  3,  the  electron  and  ion  number  densities 
will  be  relatively  low  compared  to  low  Earth  orbit.  In  this  region  the  densities  are  on  the  order  of 
10  /m  ,  and  the  plasma  is  fully  ionized  and  collisions  are  negligible.  This  means  that  all  terms 
involving  collisions  are  zero,  and  that  shearing  and  viscous  forces  are  non-existent,  so  that  the 
pressure  tensor  is  isotropic,  given  by  the  ideal  gas  law,  and  there  is  no  diffusion.  Using  these 
simplifications,  the  fluid  equations  can  be  reduced  to  the  following  equations  which  are  written 
in  conservative  form: 
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P  =  (T  -l)|ff-imT1u2J  =  pressure 
y  =  degrees  of  freedom 

Maxwell’s  equations  are  used  for  the  calculation  of  the  both  the  radiation  pattern  of  the  antenna 
and  quasi-electrostatic  effects  near  the  sheath.  Used  in  conjunction  with  the  fluid  equations,  the 
number  densities  and  velocities  are  used  to  update  the  corresponding  conduction  current  term  in 
the  curl  equation  for  the  magnetic  field. 

Maxwell  s  equations  are  solved  using  a  well-established  FDTD  (Finite  Difference  Time 
Domain)  method.  The  fluid  portion  of  the  code  is  more  complicated  since  it  incorporates 
nonlinearities  which  may  give  rise  to  a  number  of  fluid  phenomena  such  as  shocks  that  would 
render  many  finite  difference  schemes  unstable.  The  fluid  portion  of  the  code  is  solved  using  the 
method  outlined  in  [Kurganov  and  Tadmor,  2000].  This  method  involves  incorporating  a  small 
amount  of  artificial  viscosity  into  the  code  to  smooth  out  any  numerical  artifacts  that  come  about 
due  to  the  nonlinear  terms  present  in  the  fluxes.  The  artificial  viscosity  is  proportional  to  the 
maximum  characteristic  speeds  in  the  plasma  which  can  be  found  by  taking  a  Jacobian  of  the 
flux  matrices.  This  technique  has  been  proven  very  accurate  for  systems  of  conservation  laws 
[Kurganov  and  Tadmor,  2000]. 

As  in  all  finite  difference  methods,  a  mesh  is  used  to  discretize  continuous  space,  and  all 
the  derivatives  are  solved  using  samples  of  this  space.  The  configuration  for  this  model  is  shown 
in  Figure  2,  representing  a  dipole  in  three-dimensional  space.  All  fluid  variables  are  located  at 
the  cell  centers.  Electric  fields  are  located  at  the  cell  edges,  and  magnetic  fields  as  well  as  fluxes 
are  located  at  the  cell  faces. 
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As  a  result  of  the  FDTD  method  used  to  solve  for  the  radiated  fields,  E  +  H  are  staggered 
in  both  time  and  space.  The  fluid  variables  must  therefore  take  this  into  account. 

Since  the  fluid  variables  are  solved  using  a  forward  time  centered  space  (FTCS)  representation, 
the  fluid  variables  must  be  co-located  in  time  with  either  E  or  H.  We  chose  the  fluid  variables  to 
be  co-located  temporally  with  the  magnetic  field. 

Because  of  this,  both  the  old  and  new  time  steps  for  the  electric  field  must  be  stored  in  different 
arrays.  The  temporal  arrangement  for  the  different  variables  is  shown  in  Figure  3. 


Figure  2,  Mesh  and  Cell  Configuration 
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En_1  En  En+1 


Since  the  fluid  variables  U  are  solved  with  a  forward  time  scheme,  U  must  be  updated 
using  values  located  at  the  previous  time  step  and  not  halfway  in  between  like  FDTD. 
Therefore  to  solve  for  Un+,/’,  the  variables  Un  Va,Hn  Vi,  and  the  average  of  En  and  E"'1  must  be 
used.  Hence  the  need  to  store  the  old  values  of  Electric  Field. 

Figure  3.  Orientation  of  Variables  along  time  axis 


In  an  FDTD  method,  the  antenna  is  excited  by  applying  a  potential  difference  between  the  two 
conductors  of  a  dipole  as  shown  in  Figure  4.  Once  the  antenna  is  excited  (via  electric  field 
enforcement)  the  fluid  and  magnetic  field  values  are  determined  using  the  newly  acquired 
electric  field  value*  and  finally  the  updated  electric  field  is  determined  using  the  fluid  and 
magnetic  field  variables  that  were  just  calculated.  We  use  the  Fluid  code  to  find  the  antenna 
current  distribution  in  the  following  way.  The  sinusoidally  varying  potential  is  applied 
in  the  gap  between  the  conductors.  The  voltage  and  current  waveforms  propagate  along  the 
antenna  setting  up  the  correct  distributions  for  a  radiating  element.  The  electromagnetic  and 
induction  fields  are  mapped  from  the  vicinity  of  the  antenna  out  to  the  sheath  boundary  assuming 
a  1/r  variation.  The  effects  of  the  sheath  are  introduced  through  a  boundary  condition  on  the 
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1. )  Antenna  Excitation 

2. )  Solve  for  Fluid  Variables  (Number  density,  Momentum,  and  Energy) 

PR 

4. )  V  x  E  =  — —  (Solve  for  H) 

ct 

fST~\ 

5. )  Vx//  =  J  +  —  (Solve  for  E) 

6. )  Repeat  steps  1  -  5  until  convergence 


Figure  4.  Antenna  Excitation 
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radial  component  of  the  quasi-static  electric  field,  at  the  sheath  edge. 

A.2  Boundary  Conditions 

The  correct  boundary  conditions  are  perhaps  the  most  important  part  of  any  simulation. 

Boundary  conditions  in  this  model  are  made  for  the  truncation  of  both  the  FDTD  and  fluid  lattice 
as  well  as  the  boundary  conditions  for  fluxes  onto  the  antenna  surface.  Since  it  is  important  that 
both  electrostatic  and  electromagnetic  waves  do  not  reflect  back  from  the  boundary,  absorbing 
boundary  conditions  have  been  developed  for  both  the  fluid  and  Maxwell’s  equations. 

Absorbing  boundary  conditions  for  the  solution  of  Maxwell’s  equations  are  well  known. 

Perhaps  the  most  notable  and  promising  is  the  PML  (Perfectly  Matched  Layer)  formulation 
introduced  in  [Berenger,  1994].  This  requires  a  thin  layer  of  cells  around  the  FDTD  space  which 
simulate  a  highly  absorptive  media.  A  better  implementation  of  this  method  is  discussed  in 
[Roden  and  Gedney,  2000]  which  allows  for  the  absorption  of  evanescent  modes  and  is 
completely  independent  of  the  host  medium.  This  is  the  method  of  choice  for  the  EM-fluid 
model  presented  here.  Since  the  method  utilizes  a  convolution  operator,  the  equations  must  be 
approximately  linear  for  the  method  to  work.  It  has  been  shown  that  this  methods  has  an 
effective  absorption  rate  of  up  to  50-100  dB  of  the  incident  wave. 

The  PML  method  works  well  for  Maxwell’s  equations,  however  the  nonlinear  fluid 
equations  are  much  more  of  a  challenge.  Some  assumptions  must  be  made.  Equivalent  methods 
for  fluid  dynamics  are  not  yet  available  for  non-linear  equations  [Hesthaven.  1997].  However,  a 
number  of  aspects  particular  to  the  problem  at  hand  allows  for  the  use  of  this  convolutional  PML 
derivation.  Outside  of  the  sheath  region,  the  bulk  velocity  of  the  plasma  is  essentially  zero. 
Meaning  that  there  is  no  mean  fluid  flow.  All  other  quantities  are  small  relative  to  the  ambient 
values.  This  means  that  the  convective  terms  in  the  fluid  equations  are  essentially  zero,  allowing 
for  the  linearization  of  the  fluid  equations.  Boundary  conditions  for  a  zero  mean  flow  essentially 
reduce  to  those  involving  wave  equations  such  as  Maxwell’s,  thus  we  can  apply  the  same  PML 
scheme  to  the  fluid  equations.  Absorbing  boundary  conditions  must  be  implemented  for  each 
model  region  (refer  to  Figure  1)  in  order  to  prevent  contamination  from  spreading  into  each 
successive  calculation. 

Boundary  conditions  for  the  conductor  must  also  be  taken  into  account.  In  many  papers, 
it  is  assumed  that  the  particles  hit  the  walls  of  a  conductor  with  the  thermal  velocity.  The 
unidirectional  continuity  flux  of  particles  at  the  wall  is  thus  given  by: 

r  1 
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[Lymberopoulos,  Economou.  1995],  Some  papers  also  include  a  secondary  emission  coefficient 
[Meezan,  Capelli.  2001],  In  the  magnetosphere  at  L=2  in  the  magnetic  equatorial  plane,  the 
plasma  temperature  is  only  about  l/5eV,  Therefore  any  potential  in  the  kV  range  will  accelerate 
electrons  well  above  the  thermal  velocities  and  this  must  be  taken  into  account.  Therefore, 
boundary  conditions  similar  to  those  found  in  [Lymberopoulos  and  Economou,  1995 ]  are  used 
in  this  model  including  secondary  electron  emission  with  the  exception  that  the  velocity  of  the 
particles  is  actually  the  sum  of  the  thermal  velocities  into  the  wall  and  the  bulk  fluid  velocity  at 
the  wall  as  shown  in  Figure  5.  Boundary  conditions  involving  thermal  velocities  are 
implemented  for  the  continuity  and  energy  equations.  It  is  assumed  the  there  is  no  gradient  in 
any  of  the  conserved  quantities  at  the  wall  of  the  conductor. 


^  wall  ^  th  ^  bulk 

Figure  5.  Fluid  Boundary  Condition 


A.3  Antenna  Pattern 

Once  the  number  densities  reach  a  steady  state  distribution,  the  antenna  currents  and  the  currents 
within  the  sheath  region  are  used  to  determine  the  antenna  pattern  using  a  far-field 
transformation  (Fig.  6).  The  far  field  transformation  involves  the  use  of  the  radiated 
electromagnetic  fields  E  and  H,  Green’s  Theorem  and  the  Surface  Equivalence  Theorem.  The 
radiated  fields  within  the  problem  space  can  be  weighted  by  the  free  space  Green’s  function  and 
integrated  around  a  closed  surface  to  provide  the  far-field  pattern  of  the  antenna. 


9 


Imaginary  Surface  S 


j  M  t  = -fix  E 
\ 

|  J,  =nxH 


1.)  The  surface  equivalence 
theorem  is  used  to  determine 
the  equivalent  electric  and 
magnetic  fields  along  a  surface 
S  (and  hence  the  virtual 
currents  M  and  J)  surround  in; 
the  antenna. 


Far-fieid  observation  point 

f  =  {xx+$) 


2.)  Then  by  utilizing 
Green’s  functions,  and 
integrating  around  the 
surface,  the  pattern  at  a 
particular  location  can 
be  found  via 
contributions  from  each 
of  the  surface  elements. 


M  ear-field  source  point 
r-Vx+f?) 


Figure  6.  Near-field  to  Far-field  Transformation 
A.4  Requirements  of  the  full  EM-fluid  formulation 

There  are  a  number  of  requirements,  which  must  be  met  in  order  to  properly  use  the  full  EM- 
fluid  model.  In  order  for  the  boundary  conditions  to  work  effectively,  the  truncation  of  the  space 
must  occur  in  a  region  where  the  assumptions  in  which  the  boundary  equations  were  formulated 
are  satisfied.  For  instance,  one  must  be  in  the  far-field  for  the  PML  condition  to  work  for 
Maxwell’s  equations,  and  for  the  fluid  equations,  the  space  must  be  far  enough  beyond  the 
sheath  region  for  the  eold-plasma  approximation  to  be  satisfied,  i.e.  any  non-linear  equations 
must  be  well  approximated  by  linear  versions.  Since  both  Maxwell’s  and  the  fluid  equations  are 
solved  simultaneously,  the  far-field  condition  for  Maxwell’s  equation  (being  the  greater  of  the 
two)  takes  precedence. 


10 


An  option  exists  to  reduce  the  number  of  cells  necessary  for  the  simulations,  known  as  the 
AMR  (Adaptive  Mesh  Refinement)  technique.  It  involves  applying  resolution  only  where  it  is 
necessary,  such  as  around  the  radiating  structure,  where  electrostatic  effects  such  as  sheath 
development  must  be  considered.  This  method  is  useful  in  any  model  in  which  multiple 
wavelengths  exist  at  once,  such  as  the  electrostatic  sheath  model,  or  in  an  anisotropic  FDTD 
model  where  the  refractive  index  changes  with  both  frequency  and  angle  with  respect  to  the 
magnetic  field.  A  diagram  of  this  concept  is  displayed  in  Figure  7  where  the  dipole  antenna  is 
shown  in  the  middle  of  the  space.  This  method  has  been  developed  for  variables  that  are  located 
at  cell  centers  or  cell  faces  such  as  the  fluxes  and  fluid  variables  in  the  fluid  equations.  However, 
in  the  case  of  FDTD  for  Maxwell’s  equations,  the  electric  fields  are  located  at  the  cell  edges,  to 
allow  for  an  accurate  curl  representation.  The  current  AMR  method  is  not  applicable  to  this 
case,  however  we  estimate  that  appropriate  modifications  to  the  AMR  method  for  use  in  our 
simulations  could  be  completed  in  3-6  months. 

Concept  of  Adaptive  Mesh  Refinement  (AMR)  for  Antenna  Problem 
Largest  Cell  Size  ~  100m 
Smallest  Cell  Size  ~  ,01m 

This  means  that  the  problem  requires  about  14  levels  of  refinement  if  each  level  of 
refinement  is  considered  to  be  a  reduction  by  2  in  cell  size 


Figure  7.  AMR  (Adaptive  Mesh  Refinement) 
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A.5  Poisson  Model  for  the  Sheath 

.  At  potentials  in  the  kV  range,  the  sheath  may  be  on  the  order  of  a  few  meters  in  radius.  Since 
the  driving  frequency  of  the  antenna  will  generally  be  much  less  than  both  the  electron  and  ion 
plasma  frequencies,  the  sheath  should  reach  a  quasi-steady  state  at  each  point  in  the  sinusoidal 
cycle.  Therefore  at  each  point  in  the  sinusoidal  variation,  the  sheath  should  expand  and  collapse 
with  the  potential  of  the  driving  elements.  The  characteristics  of  the  sheath  dimensions  and  the 
dynamics  of  the  charged  particles  within  the  sheath  can  then  be  determined  through  a  simulation 
involving  a  Poisson  model. 

The  Poisson  model  involves  solving  the  fluid  equations  with  electric  field  source  terms 
determined  through  the  use  of  Poisson’s  equation. 

V20>  =  — 

*0 


Poisson’s  equation  is  solved  using  a  fast-sine  transform,  which  falls  under  the  category  of 
rapid  elliptic  solvers.  This  technique  proved  to  be  the  most  efficient  method  for  solving  elliptic 
equations  of  this  sort.  To  place  conducting  objects  at  specified  potentials,  it  is  necessary  to 
incorporate  a  capacity  matrix  method,  which  allows  sources  internal  to  the  system  to  be  placed 
inside  the  space  and  provide  a  correct  solution  of  Poisson’s  equation  [Hockney  and  Eastwood, 
1981]. 

In  using  a  Poisson  type  technique  it  is  assumed  that  the  potential  is  constant  along  each 
dipole  element.  This  is  equivalent  to  assuming  a  triangular  current  distribution  along  the 
antenna,  which  is  the  distribution  expected  for  antennas  whose  length  is  much  shorter  than  a 
wavelength.  In  the  Poisson  approach,  the  potentials  are  applied  on  each  of  the  conductors 
independently  as  shown  in  Figure  8,  however  the  magnitude  of  the  applied  potentials  is  assumed 
to  vary  sinusoidally. 
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Two  constant  potentials  of +V  and  - 
V  are  applied  across  the  entire 
conductor  in  this  electrostatic  case 
(The  voltage  waveform  is  a  constant 
across  the  length  of  the  conductor) 


Figure  8.  Poisson  Model 

The  timestep  constraint  for  Poisson  type  models  is  well  documented  as  being  the  minimum  of 
either  the  inverse  of  the  plasma  frequency,  that  which  is  related  to  the  maximum  characteristic 
speed  in  the  plasma  based  on  thermal  effects  and  bulk  velocities,  or  the  dielectric  relaxation  time 
if  the  plasma  is  eollisional  which  is  often  the  most  restrictive  of  the  timesteps.  Since  the  plasma 
in  the  current  formulation  is  collisionless,  the  last  constraint  does  not  apply.  For  the  operating 
conditions  of  this  model,  the  maximum  characteristic  speeds  in  the  plasma  are  generally  the 
determining  factor.  In  our  case,  the  timestep  corresponding  to  the  thermal  velocities  are  more 
restrictive  than  those  having  to  do  with  plasma  oscillations.  And  for  a  0.2  eV  plasma,  and 
applied  potentials  on  the  antenna  of  1 000Y  or  more,  the  bulk  velocities  almost  always  exceed  the 
thermal  velocities  in  a  collisionless  regime. 

Trying  to  resolve  a  steady  state  solution  on  small  grids  with  very  low  frequency  waves 
poses  a  challenging  problem  at  the  present  time,  but  the  complexity  can  be  reduced  to 
manageable  levels  using  the  AMR  technique.  An  alternative  approach  exists  which  involves  the 
solution  of  Poisson’s  equation  for  a  slowly  varying  applied  potential  over  subsections  of  the 
large  antenna  whose  lengths  are  much  larger  than  the  sheath  radius.  By  slowly  varying  we  mean 
that  the  frequency  of  the  applied  potential  is  much  smaller  than  either  the  electron  or  ion  plasma 
frequency,  which  will  generally  be  the  case  for  the  proposed  VLF  radiating  system.  In  this  way 
the  sheath  characteristics  as  a  function  of  applied  potential  could  be  calculated  along  the  entire 
antenna.  Using  these  sheath  characteristics  developed  in  the  Poisson  model,  it  would  then  be 
possible  to  implement  a  dynamic  sheath  into  an  FDTD  approach.  Essentially,  this  involves 
merging  the  electrostatic  and  electromagnetic  models  as  outlined  in  Figure  1 . 
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A.6  Anisotropic  Cold  Plasma  FDTD  Model 


In  region  3  of  Figure  1,  the  simulation  model  involves  the  solution  of  Maxwell’s 
equations  in  a  cold  plasma  and  the  propagation  of  ELF/VLF  waves  in  this  inhomogeneous, 
anisotropic  medium.  Once  the  currents  in  regions  1  and  2  of  Figure  1  are  specified,  the  cold 
plasma  model  can  then  be  used  to  determine  the  radiation  fields  produced  by  these  currents.  This 
is  an  alternative  method  to  that  shown  in  Figure  6.  The  model  involved  utilizes  the  method 
outlined  in  [Lee  and  Kalluri,  1999],  This  includes  the  solution  of  the  following  equations: 


Vx£ 


-Mo 


dH 


dt 
BE 


VxH  =  e0—+J 
0  dt 

dJ 

—  +  vJ  =  e0a)p  (r,t)E  +  (ob{r,t) x  J 


Here,  s0  is  the  free  space  permittivity,  p0  is  the  free  space  permeability,  v  is  the  collision 
frequency,  o>p  is  the  plasma  frequency,  and  coD  is  the  electron  gyro  frequency.  Warm  plasma 
effects  can  be  introduced  through  a  pressure  term  when  necessary  [Young  and  Brueckner,  1994], 
These  equations  can  be  used  to  accurately  model  wave  propagation  outside  of  the  sheath  and  hot 
plasma  regions.  To  demonstrate  the  validity  of  the  model,  we  present  a  series  of  plots  showing  a 
simple  cold  plasma  regime  with  operating  points  specified  on  the  w-k  diagram  shown  in  figure  8. 


.  0)  vsK 

X10 


Figured  o-k diagram 
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The  model  however,  was  thoroughly  tested  in  1-d  for  points  on  the  ordinary  and  parallel 
branches  as  well.  The  model  accurately  depicted  cutoff  regions  and  displayed  appropriate  wave 
numbers  for  the  test  cases. 

In  the  following  section,  the  validity  of  the  multidimensional  anisotropic  model  in  a  cold 
plasma  is  demonstrated.  Both  free  space  propagation  and  that  within  a  plasma  are  presented.  A 
2-dimensional  coordinate  system  is  set  up  to  allow  for  calculation  of  both  the  free  space  and 
extraordinary  modes  on  the  w-k  diagram.  In  the  following  2-d  models,  the  antenna  is  a  Vik 
dipole  with  Hz,  Ex,  and  Ey  as  the  field  configuration  and  a  static  magnetic  field  directed  in  the  z- 
direction.  The  antenna  elements  are  in  the  x-y  plane  with  the  antenna  oriented  along  the  y- 
direction  and  infinite  in  extent  in  the  z-direction  as  shown  in  Figure  10.  The  antennas  are 
excited  via  a  potential  difference  across  the  gap  between  the  elements  in  the  y-direction  (Ey 
excitation). 

A  plot  of  the  Fourier  transform  of  the  spatial  frequency  is  also  displayed.  This  will 
provide  a  comparison  for  the  cold  plasma  case.  The  following  figures  will  show  the  same  dipole 
configuration  for  an  antenna  operating  in  the  plasma  at  the  location  specified  by  the  operating 
point  in  Figure  9.  The  2-d  simulation  plots  shown  with  the  exception  of  the  free  space  diagram 
will  be  those  representing  the  extraordinary  branch  of  the  w-k  diagram. 


Figure  10.  2-d  antenna  simulation 

Figure  1 1  shows  the  plot  of  the  free  space  Ey  field  after  1 7  cycles  at  491kHz.  The  spatial 
Fourier  transform  at  this  time  step  is  plotted  in  Figure  12  and  a  zoom  in  of  the  FFT  is  displayed 


Spatial  frequency  1  /meters  meters 


in  Figure  13.  At  an  operating  frequency  of  491kHz  in  free  space,  the  wavelength  should  be 
m  and  the  spatial  frequency  should  be  about  .0016  1/m. 


Figure  11.  Ey  field  in  free  space 
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Figure  12.  FFT  of  Ev  field 
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Figure  13.  Enlargement  of  Ey  field  FFT 

It  can  bee  seen  from  the  plot  that  the  spatial  frequency  is  indeed  1 .6e-3  1/m  as  one  would  expect. 
Figure  14  shows  the  propagation  of  the  Ey  field  for  the  free  space  case  over  a  period  of  1  cycle 
of  the  operating  frequency. 

To  test  the  validity  of  the  anisotropic  model,  the  extraordinary  branch  at  a  frequency  of 
625  kHz  was  modeled.  At  this  frequency,  the  refractive  index  is  equal  to  2.  Therefore  the 
wavelength  should  be  !4  of  the  free  space  value.  For  the  anisotropic  case,  this  spatial  frequency 
should  be  .0042  1/meters.  The  same  set  of  plots  is  made  for  the  Ey  field  in  the  cold  plasma. 
Figure  15  shows  the  Ey  field.  Figure  16  and  17  show  the  FFT  and  enlargement  of  the  FFT  for 
this  field  after  17  cycles.  Once  again  good  agreement  is  obtained.  Figure  18  is  another  sequence 
of  the  Ey  field  over  1  period  for  the  anisotropic  case. 
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Figure  17.  Enlargement  of  Ey  field  FFT  in  Cold  Plasma 
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The  final  ease  is  a  three  dimensional  run  made  at  the  location  defined  as  the  operating 
point  on  Figure  19  which  corresponds  to  15kHz.  In  this  case  all  components  of  E,  H,  and  J  are 
used  in  the  simulation.  The  model  is  being  run  at  a  frequency  corresponding  to  the  whistler 
mode.  The  same  configuration  of  antenna  and  static  magnetic  field  as  in  the  2-d  case  are  used. 
The  static  magnetic  field  is  oriented  in  the  z-direction  and  the  antenna  along  the  y-direction  in 
the  x-y  plane.  The  only  difference  is  that  the  z-dependence  of  the  variables  is  added.  It  can  be 
seen,  that  as  expected,  the  fields  are  directed  toward  the  static  magnetic  field,  which  is  in  the  z 
direction.  Sequences  of  slides  for  the  three  slices  through  the  antenna  are  shown  on  the  following 
pages  for  the  Ey  component  which  is  the  field  of  excitation  on  the  antenna.  Evanescent  modes 
are  seen  in  the  x-y  plane,  and  propagation  is  seen  within  a  very  small  cone  centered  along  the  z- 
direction  as  expected.  Since  in  FDTD  it  is  necessary  to  resolve  the  smallest  wavelength  in  the 
system,  collisions  were  added  to  control  the  effects  at  the  resonance  cone  where  X— >  0. 


Figure  19.  <a-k  diagram 


Figure  21.  Ey  field  in  X-Y  Plane  (Evanescent) 


Figure  22.  Ey  field  in  X-Z  Plane.  Propagation  along  Z. 
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B  Antenna  Input  Impedance. 

A  second  main  Stanford  goal  has  been  to  determine  the  first  order  effects  of  the  dipole  antenna 
plasma  sheath  on  the  antenna  input  impedance  and  to  devise  methods  to  control  these  effects.  In 
this  study  we  have  used  numerical  solutions  to  the  combined  Poisson’s  and  Lorentz’s  equations 
to  determine  the  sheath  characteristics  in  two  dimensions.  We  treat  each  antenna  element  as  an 
infinitely  long  cylinder  and  determine  the  sheath  capacitance  per  meter  of  antenna  length.  Our 
approach  is  similar  to  that  of  LaFramboise  and  Sonmor  [1993] 

Working  from  our  2-D  model,  we  have  discovered  a  new  method  of  controlling  both  the  sheath 
dimensions  and  the  sheath  capacitance  when  the  antenna  is  used  to  radiate  VLF  waves  in  the 
magnetospheric  plasma.  This  is  an  important  topic  since  the  sheath  capacitance  is  responsible  for 
a  large  portion  of  the  antenna  input  impedance. 

To  illustrate  the  new  control  method,  we  first  show  in  Figure  26  the  predictions  of  our  2- 
D  model  for  the  sheath  radius  r  as  a  function  of  applied  antenna  potential.  It  can  be  seen  that  the 
sheath  radius  varies  from  a  few  centimeters  for  no  applied  voltage  up  to  a  few  meters  at  1000 
Volts.  This  is  roughly  the  range  of  voltages  that  would  be  applied  to  the  dipole  antenna  during 
normal  transmitting  operations. 

Since  the  sheath  capacitance  is  a  function  of  r,  we  can  expect  similar  variations  in  the 
sheath  capacitance  as  the  applied  voltage  is  varied.  Figure  27  shows  the  predicted  sheath 
capacitance  as  a  function  of  r  for  a  long  dipole  antenna.  It  varies  widely  from  approximately  25 
pF/m  at  low  applied  voltage  to  approximately  5  pF/m  in  the  kV  range. 

For  a  long  dipole  antenna,  the  sheath  capacitance  is  responsible  for  a  major  part  of  the 
input  impedance  of  the  antenna.  However  if  we  drive  the  antenna  with  a  sinusoidal  voltage  of 
roughly  1000  V  it  is  clear  that  the  sheath  capacitance  will  vary  greatly  during  an  RF  cycle.  Thus 
it  is  not  possible  to  define  an  input  impedance  in  a  normal  sense  for  the  antenna.  This  is  a  major 
problem.  If  the  antenna  cannot  be  tuned,  the  radiation  efficiency  of  the  spacecraft  VLF 
transmitting  system  will  be  very  low  and  mission  success  will  be  compromised. 

Using  our  2-D  sheath  model,  we  have  discovered  that  the  problems  inherent  in  a  widely 
varying  sheath  capacitance  can  be  mitigated  by  applying  a  DC  bias  to  both  elements  of  the  dipole 
antenna.  The  bias  can  be  either  positive  or  negative.  A  large  positive  bias  voltage  can  be  applied 
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using  a  power  source  and  an  electron  gun  on  the  spacecraft.  Negative  current  carried  to  the 
antenna  by  electrons  is  collected  and  expelled  from  the  spacecraft  using  the  power  source 
and  an  electron  gun.  The  maximum  amount  of  bias  achievable  through  this  method,  or  similar 
methods,  will  be  a  subject  for  future  study. 

To  illustrate  the  concept,  we  assume  that  a  positive  DC  voltage  bias  of  5000  V  is  applied 
to  both  halves  of  a  dipole  antenna  and  the  AC  voltage  is  then  swept  from  -2000  V  to  +  2000  V 
about  the  bias  voltage.  Results  from  our  2-D  sheath  model  are  shown  in  Figure  28.  The  Figure 
shows  the  sheath  capacitance  as  a  function  of  the  sweeping  voltage.  It  can  be  seen  that  the  sheath 
capacitance  varies  only  by  approximately  10%  over  the  4000  V  sweep  range.  These  results 
suggest  that  this  biasing  technique  can  be  very  effective  in  controlling  the  dipole  sheath 
capacitance. 

For  a  positive  bias  of  5000  V,  our  2-D  model  predicts  that  the  electron  current  to  the 
antenna  would  be  approximately  20  microAmps/m.  Thus  a  100  m  tip-to-tip  dipole  would  draw 
about  2  mA,  and  the  electron  gun  system  would  expend  about  10  W  in  removing  the  arriving 
negative  charge  from  the  antenna.  In  addition  there  would  be  some  power  necessary  to  run  the 
electron  gun.  The  use  of  an  ion  gun  along  with  a  large  negative  bias  might  also  be  effective.  The 
ion  current  arriving  at  the  negatively  biased  antenna  would  be  much  smaller  than  the  electron 
current  expected  for  a  positive  bias  of  the  same  magnitude.  Thus  the  power  expended  by  the  ion 
gun  system  in  removing  the  arriving  positive  charge  would  be  smaller,  but  more  energy  might  be 
required  to  expel  the  ions  from  the  ion  gun  because  of  their  larger  mass.  There  are  clearly  a 
number  of  questions  that  need  to  be  addressed  concerning  this  biasing  technique.  Nevertheless 
this  method  appears  very  promising. 

In  the  coming  three  year  period  we  propose  to  determine  the  feasibility  of  using  high 
power  electron  or  ion  guns  on  the  transmitting  spacecraft  to  stabilize  the  sheath  capacitance 
and  allow  straightforward  tuning  of  the  dipole  antenna.  We  will  also  examine  the  possibility  of 
using  an  active  tuning  network  that  compensates  in  real  time  for  the  sheath  capacitance  variation 
during  each  RF  cycle. 


100 


200  300 

Sheath  Radius  (cm) 


400 


50 


Figure  24.  Antenna  voltage  and  resultant  plasma 
sheath  radius  for  two  values  of  plasma  temperature. 
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Figure  25,  Sheath  capacitance  as  a  function  of  sheath  radius. 
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Antenna  Voltage  (V) 

Figure  26.  Sheath  capacitance  as  a  function  of  antenna 
voltage  for  a  positive  bias  voltage  of  5000  V. 
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C.  Other  Accomplishments 

Additional  goals  of  the  Stanford  effort  have  been:  1)  the  determination  of  nominal  antenna 
requirements  from  wave-particle  interaction  models  that  are  required  to  maximize  particle  decay 
rates  for  targeted  populations,  2)  the  determination  of  the  effects  of  Landau  damping  on  the 
waves  radiated  from  the  antenna,  3)  the  determination  of  the  energetic  electron  precipitation  due 
to  magnetospherically  reflected  whistlers,  4)  the  development  .of  a  global  overview  of  the  total 
radiated  VLF  power  from  space  necessary  to  match  the  effects  of  presently  operating  ground 
based  VLF  transmitters,  and  5)  the  determination  of  the  VLF  input  impedance  of  .large  loop 
antennas  in  the  magnetosphere.  All  of  these  topics  must  be  considered  when  evaluating  the 
feasibility  of  radiation  belt  control  using  space  based  VLF  transmitting  systems. 

C.l .  Antenna  Length 

For  a  variety  of  spacecraft  orbits,  we  have  determined  the  optimum  dipole  antenna  length  for 
maximizing  the  radiation  of  ELF/VLF  waves  with  the  correct  wave  normal  angles  for 
precipitating  target  electrons  in  the  1-10  Mev  energy  range.  This  work  is  reported  in  [Inan  et  al. 
,2002}. 

C  2.  Landau  Damping 

We  have  used  data  from  the  POLAR  spacecraft  to  determine  the  characteristics  of  energetic 
electrons  in  the  100  eV  -  2  keV  energy  range  within  the  radiation  belts.  These  characteristics 
were  then  used  to  assess  the  Landau  damping  that  would  affect  the  amplitude  of  the  radiated 
ELF/VLF  waves  as  they  propagated  through  the  belts.  Knowledge  of  this  damping  and  its 
variation  with  wave  normal  angle  is  critical  in  the  evaluation  of  the  feasibility  of  the  mission. 
This  work  is  reported  in  [Bell,  et  al,  2002.] 

C.3.  Particle  Precipitation  by  Magnetospherically  Reflected  Whistlers 

Using  Landau  damping  rates  from  the  foregoing  study,  we  have  determined  the  energetic 
electron  precipitation  that  is  produced  by  magnetospherically  reflected  whistlers.  This 
precipitation  mode  is  one  which  we  can  emulate  using  spacebased  ELF/VLF 
transmitters. [Bortnik  et  al.,  2002], 
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C.4.  Global  Overview 

We  worked  with  Tony  Tether  of  ARPA  to  provide  a  global  overview  of  the  number, 
placement ,  and  radiated  power  of  a  fleet  of  ELF/VLF  radiating  spacecraft  which  would  be 
required  to  precipitate  energetic  electron  fluxes  equal  to  those  presently  being  precipitated  by 
ground  based  VLF  transmitters.  This  work  is  reported  in  [I nan  et  ah,  2002] 

C. 5.  The  VLF  Input  Impedance  of  Large  Loop  Antennas 

We  determined  the  ELF/VLF  radiation  characteristics  of  large  (100  m  radius)  loop 
antennas  deployed  from  spacecraft  in  both  low  and  high  altitude  orbits.  In  general,  loops  used  for 
radiating  VLF  waves  must  be  roughly  of  this  size  in  order  to  radiate  sufficient  power.  This  type 
of  antenna  would  appear  to  have  an  important  advantage  over  dipole  antennas  in  low  Earth  orbit, 
since  a  loop  antenna  with  constant  current  does  not  develop  a  large  plasma  sheath  that  would 
affect  its  input  impedance.  However  our  results  suggest  that  even  small  spatial  variations  of 
current  around  the  loop  can  greatly  affect  the  input  impedance.  This  work  is  reported  in  [Bell  et 
ah,  2002], 

D.  Joint  Work 

During  the  duration  of  his  Grant  we  have  worked  jointly  on  various  topics  with 
AFRL/VSB  personnel.  Our  work  with  Dr.  J.  Albert  has  led  to  the  presentation  of  a  joint  paper 
at  the  Fall  2001  AGU  meeting  on:  “Optimal  VLF  wave  parameters  for  pitch  angle  scattering  of 
trapped  electrons”  and  to  the  submission  of  a  joint  paper  to  the  Journal  of  Geophysical  Research 
concerning  controlled  precipitation  of  radiation  belt  particles  [  Inan,  et  al,,  2002]. 
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